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Abstract. In recent years, a number of finite element methods have been formulated for the 
solution of partial differential equations on complex geometries based on non-matching or overlapping 
meshes. Examples of such methods include the fictitious domain method, the extended finite clement 
method, and Nitsche's method. In all of these methods, integrals must be computed over cut cells or 
subsimplices which is challenging to implement, especially in three space dimensions. In this note, 
we address the main challenges of such an implementation and demonstrate good performance of a 
fully general code for automatic detection of mesh intersections and integration over cut cells and 
subsimplices. As a canonical example of an overlapping mesh method, we consider Nitsche's method 
which we apply to Poisson's equation and a linear elastic problem. 
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1. Introduction. A fundamental problem in computational science is to develop 
methods for the solution of partial differential equations on domains containing one or 
several objects that may have complex or time-dependent geometry. One approach to 
attacking this problem is to allow overlapping meshes where a mesh representing an 
object is allowed to overlap a background mesh representing the surroundings of the 
object; see for instance Chesshire and Henshaw [22], Mayer et al. [50], Yu [60], Zhang 
et al. [61] for various applications. 

In solid mechanics, overlapping meshes may be used to represent materials con- 
sisting of elastic objects inserted into a surrounding elastic material of another type 
[37], and in fluid-structure interaction, an overlapping mesh may be used to represent 
an elastic body immersed in a fluid represented by a fixed background mesh [7, 49, 60]. 
Another common application [22, 25, 26] is found in mesh generation where a com- 
plicated geometry such as, for example, a pipe junction, may be decomposed into 
simpler parts and one unstructured tetrahedral mesh is created for each part. These 
components may then be stored, reused and recombined in applications by using an 
overlapping mesh technique. 

Overlapping mesh techniques are of particular interest in simulations that involve 
moving objects. For such problems, overlapping mesh techniques are an attractive 
alternative to ALE techniques. The main advantage is that by using an overlapping 
mesh technique, one avoids deformation of the mesh that may lead to deterioration of 
the mesh quality and ultimately force remcshing. This is of particular importance in 
the simulation of fluid-structure interaction where the topology of the fluid domain 
may change due to deformation of the solid. 

The main focus of this work is on the general algorithms and efficient implemen- 
tation that is required to handle complex problems posed on overlapping meshes in 
three space dimensions. Many of the presented algorithms and the tools developed are 
of interest and use for the implementation of various overlapping mesh techniques. 
To make the discussion concrete, we here focus mainly on Nitsche's method. In 
Hansbo et al. [33], a consistent finite element method for overlapping meshes based 
on Nitsche's method was introduced and analyzed. The basic idea is to construct 
a finite element space by taking the direct sum of the space of continuous piece- 
wise polynomial functions on the overlapping mesh and the restriction of the space 
of continuous piecewise polynomial functions to the complement of the overlapping 



mesh, and then impose the interface conditions using Nitsche's method. It was shown 
that this approach leads to a stable method of optimal order for arbitrary degree 
polynomial approximation. 

The main challenge in the implementation is to compute the intersection between 
the overlapping and the background mesh. The result is a set of cut cells which may 
be arbitrarily complex polyhedra. These arise as the result of subtracting from the 
tetrahedra of the background mesh a set of tetrahedra from the overlapping mesh. 
By adopting algorithms and search structures from the field of computational geom- 
etry, we show how these issues can be handled in an efficient manner. Furthermore, 
one needs to compute integrals on the resulting polyhedra. This can be carried out 
efficiently based on an application of the divergence theorem in combination with 
potentials to represent an integral on the three-dimensional polyhedron as a sum of 
one-dimensional integrals on its edges. 

The presented algorithms and implementation are relevant for several other types 
of related methods, including the extended finite clement method (XFEM) [31, 58], 
non-fitted sharp interface methods [14, 32] , and mesh-tying techniques [25] . 

1.1. Major contributions of this paper. Our work consists of several contri- 
butions. We identify the major techniques used in the implementation of overlapping 
mesh methods and related methods. We further identify useful data structures and 
algorithms from the field of computational geometry. As part of our work, existing 
computational geometry libraries such as CGAL [1] and GTS [2] have been wrapped into 
the general purpose finite element library DOLFIN [43] and into an extension library 
on top of it, thereby making these algorithms and data structures more easily acces- 
sible to the finite clement community. Based on our implementation, we demonstrate 
for the first time a highly efficient implementation of Nitsche's method on overlap- 
ping meshes for several problems posed in three spatial dimensions, thus opening the 
possibility of employing Nitsche-based overlapping mesh methods for challenging 3D 
problems such as fluid-structure interaction or domain-bridging problems. 

1.2. Outline of this paper. In Section 2, we review Nitsche's overlapping mesh 
method for a model problem and present in Section 3 the techniques and algorithms 
we have developed for efficient implementation of Nitsche's method in three space 
dimensions. The corresponding implementation and data structures are described in 
Section 4. Finally, we present in Section 5 numerical examples that demonstrate the 
convergence of the numerical solution as well as the scaling of the work required to 
compute mesh intersections relative to standard finite element assembly on matching 
meshes. 

2. Nitsche's method on overlapping meshes. We here review Nitsche's 
method for a simple model problem posed on two overlapping meshes. 

2.1. Model problem. Let 51 = (ill U ri 2 )° be a domain in M 3 , with boundary 
<9f2o, consisting of two (open and bounded) subdomains J7i and ri 2 separated by the 
interface T = dili fl dfl2- We consider the following elliptic model problem: find 
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Fig. 2.1. Two domains f!i and Q2 separated by the common interface F = dQi D 9^2- 
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Wc here use the notation Vi — v\ci i for the restriction of a function v to the subdomain 
fli for « = 1,2, n is the unit normal to T directed from into Oi, and [v] = 
v 2 — v 1 denotes the jump in a function over the interface T. In addition, the boundary 
OUq is divided into two subdomains dflo,u and <9S!o,n where Dirichlet and Neumann 
boundary conditions are applied, respectively. 

2.2. Finite element formulation. We consider a situation where a background 
mesh 7o is given for f2 = (^1 U ^2)° and another mesh T2 is given for the overlapping 
domain 1^2 (see Figure 2.1). Both meshes are assumed to consist of shape-regular 
tetrahedra T. The mesh 71 covering Oi = (£1 \ ^2)° is constructed by 

Ti = {THQi :Te7BA|TnQi| >0}. (2.6) 

Note that the 71 consists of both regular and cut elements since the mesh To is not a 
conform tetrahedralization of the subdomain Qi in general. For a cut element T(~l fi, 
the degree of freedoms are defined to be the same as for the original uncut element 
T. This is illustrated in Figure 2.2. 

We let Vh i denote the space of continuous piecewise fixed-order polynomials on 
Ti that vanish on d^li n dtt ,D f° r i = 0, 2. The space Vh,\ is constructed as the 
restriction to Oi of functions in Vh,o; that is, Vh.\ = ■ v € Vh,o}- We may then 

define a finite element space for the whole domain fi by 

W h = V htl Q)V h ,2. (2.7) 

Nitsche' method proposed by Hansbo et al. [33] then takes the form: find Uh <G Wh 
such that 

a(u h ,v)=l{v) VweWfc, (2.8) 
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Fig. 2.2. The mesh 7i for Oi = (f!o\^2)° is constructed as the set of cells of To not completely 
overlapped by the cells of T2 ■ 



where 
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a(u,«)=V / Vu-Vvdx (2.9) 
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(Vu-n)[u]d5- / (Vvn)[u]dS + 7 / h~ l [u] ■ [v] dS, 
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Consistency Symmctrization Penalty /Stabilization 



l(v) = / fvdx+ / gvds. (2.10) 

Here, 7 is a positive penalty parameter and the average (Vv ■ n) is chosen to be 
the one-side derivative (Vt> • n) = Vt>2 • n according to Hansbo et al. [33]. But as 
pointed out by the authors, any convex combination of the normal derivatives leads 
to a consistent formulation. 

2.3. Summary of theoretical results. When the penalty parameter 7 in (2.9) 
is chosen large enough, the form a(-, •) is coercive on the discrete space Wh and one 
can derive optimal order a priori error estimates in both the energy norm and the L 2 
norm for polynomials of arbitrary degree p. The estimates take the form 

^IIV^-^H^ + ^ll^lllJ < C7^H|n , P+ i, (2.11) 

<Ch p+l \\u\\ Qo , p+1 . (2.12) 

See Hansbo et al. [33] for details. Here, || • \\n .s denotes the standard Sobolev norm 
of order s > on f^o- A posteriori error estimates and adaptive algorithms are also 
presented in Hansbo et al. [33]. 

As observed by Areias and Belytschko [6], Nitsche's method for an arbitrary 
cutting interface as in Hansbo and Hansbo [32, 34] can be reinterpreted as a particular 
instance of an extended finite element method. In contrast, the presented formulation 
on overlapping meshes lacks such a reinterpretation since two unrelated meshes are 
involved and therefore the discontinuity across the interface cannot be modeled by 
the enrichment of degrees of freedom within a single cell. Nevertheless, both methods 
share some features, especially with regard to their implementation. 



3. Techniques and algorithms. The main challenges in the implementation of 
Nitsche's method on overlapping meshes arise from the geometric computations which 
are necessary to assemble the discrete system associated with (2.8)-(2.10). Naturally, 
these geometric computations are more involved in 3D than in 2D. Similar challenges 
are encountered in related methods such as the extended finite element method [31], 
but also for the simulation of contact mechanics [59]. Different solutions have been 
proposed, see for example Sukumar et al. [58] in the case of extended finite elements 
methods, and Yang and Laursen [59] for contact problems. Here, we take another 
approach which efficiently solves the issues in the case of the Nitsche overlapping 
mesh method. In the following, we will first discuss the challenges and their remedies 
in general terms, and then return to the specific details of our implementation in 
Section 4. 

3.1. Assembly on overlapping meshes. The main implementation challenges 
arise from the fact that the interface T can cross the overlapped mesh To in an arbitrary 
manner, which has two consequences. First, the definition of the finite element space 
(2.7) involves the restriction of the function space Vh,o, defined on the full background 
mesh 7o of £Iq, to the domain The restriction results in non-standard clement 
geometries along the interface T as it can be observed from the definition 2.6 of 
7i- Secondly, the weak imposition of the interface conditions (2.2) and (2.3) by the 
interface integrals in (2.9) involves finite element spaces defined on two unrelated 
meshes. Both make the assembly challenging to implement, compared to a standard 
finite element method. 

To better understand what kinds of challenges the Nitsche overlapping mesh 
method adds, we briefly recall the general theme of the assembly routine as it is real- 
ized nowadays in many finite element frameworks; see, e.g., Bastian et al. [12], Logg 
[41]. Nitsche's method is closely related to the classical discontinuous Galerkin (DG) 
methods, which makes it natural to depart from the assembly in the DG case. A 
detailed description of finite element assembly for discontinuous Galerkin methods 
can be found in 01gaard et al. [55]. We assume a variational problem of the following 
form: find Uh G Vu such that 



where a and I are bilinear and linear forms, respectively, and Vh and Vh are the discrete 
trial and test spaces, respectively. For simplicity, we here make the assumption Vh = 
Vh- The solution of the variational problem (3.1) may then be computed by solving 
the linear system 



for U E M. N , with stiffness matrix Ajj — a(4>j,<f>i) and load vector b] — l(<fii). Here, 
{<A/}/Li denotes a basis for Vh- During assembly of the tensors A and b, one usually 
iterates over all cells T e T and computes the contributions from each cell as the cell 
tensors Afj = a(<pj , 4>J) and bf = l(4>i) where {<ff }™ =1 is the local finite element basis 
on T (the shape functions). The element tensors are then scattered into the global 
tensors A and b by adding the entries according to a given local-to-global mapping 
It '■ i i— > 1 ■ One may similarly add contributions from each facet F of the set of 
boundary facets d e T (the set of exterior facets) and, for the implementation of a DG 
method, from each facet F of the set of interior facets d{I '. Algorithm 3.1 summarizes 
the parts of the standard assembly algorithm relevant for our discussion. 



a(u h ,v) = l(v) Vv e Vh 



(3.1) 



AU = b, 



(3.2) 
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Algorithm 1 Standard finite element assembly of a bilinear form a(-, •). The local 
finite element spaces are denoted by Vh(T),Vh{F) etc. 

A = Q 

for each T e T: 

J(T) = {l,...,dim(Vfc(T))}x{l,. 
for each (i, j) e X(T) : 

for each (i, j) e I(T) : 

for each F e <9 e T: 

2 e (F) = {l,...,dim(V r fc e (F))}x{l, 
for each (i, j) e I e (F): 

Al=a{4>f 
for each (i, j) e I e {F): 

A i%(i),i° F {j)+ = A Ij 
for each F € 9,T: 

X i (F) = {l,...,dim(t^(F))}x{l, 
for each (i, j) e Ii(F): 

for each (i, j) e %i{F): 



We now consider how the standard assembly algorithm must be modified in the 
case of Nitsche's method on overlapping meshes. We first note that the tessellation 
7o of the background domain Q may be decomposed into three disjoint subsets: 

To = 7B,iU7B,2U7B,r, (3.3) 

where T ,i = {T e % : T C Hi}, T , 2 = {T e % : T C H 2 }, and T , r = {T e T : 
|T n > 0, i = 1,2} denote the sets of not, completely and partially overlapped 
cells relative to 1^2, respectively. 

Integrals over the cells of 7o,i can be assembled using a standard assembly al- 
gorithm. Furthermore, integrals over the cells of 7o,2 need not be assembled (the 
corresponding contributions will be assembled over T2 )■ However, assembly must be 
carried out over 7o,r> the partially overlapped cells of To- For the model problem 
(2.8)-(2.10), this requires the evaluation of integrals of the type 

J V<£j • Vcpf dx and J f<ff dx, (3.4) 

on cut elements (polyhedra) P = TO fl\ where T e To- Examining (2.8)-(2.10), we 
further note that we must assemble the terms 

- f • n>[&] dS- f (V& • n)[^] dS + j f ■ [&] dS. (3.5) 

This poses an additional challenge, since the integrands involve products of trial and 
test functions defined on different meshes. Furthermore, the interface T consists of a 
subset of the boundary facets d e T2 of T2, but each such facet may intersect several 
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..,dim(T4(T))} 



...,dim(^(F))} 



...,dim(^(F))} 




Fig. 3.1. The overlapping meshes in the intersection zone. The filled circles represent the 
degrees of freedom on each element (here for a piecewise linear approximation on T2 o,nd piecewise 
quadratic on To). The interface T is partitioned so that each part intersects exactly one cell Tl of 
the overlapping mesh T2 and one cell T ; ° of the background mesh To ■ 



cells of 7o- We therefore partition each facet on V into a set of polygons \Tki} such 
that each polygon Tki intersects exactly one cell T% of the overlapping mesh T2 and 
one cell X)° of the background mesh To- This is illustrated in Figure 3.1. Assembly 
may then be carried out by summing the contributions from each polygon r ki- 
ln summary, we identify the following main challenges in the implementation of 
Nitsche's method on overlapping meshes: 

1. collision detection: to determine which cells are involved in the intersection 
between the two meshes; 

2. mesh intersection: to compute the cut cells of the background mesh 7o rep- 
resented by the polyhedra {P®} and the intersection interface represented by 
the polygons {r fc/ }; 

3. integration on complex polyhedra: to compute integrals on the polyhedra {P®} 
(and the polygons {T k i}). 

In the following, we discuss these challenges in some detail and introduce concepts, 
data structures and algorithms to handle them in an efficient manner. We emphasize 
that the proposed solutions are not limited to the implementation of Nitsche-type 
methods, but may also be used for the implementation of related overlapping mesh 
methods. 

3.2. Collision detection. Topological relations between entities of a single 
mesh can be described by concepts of connectivity or mesh incidence as presented 
in Logg [42] and Bastian et al. [11]. To represent the topological relation between two 
overlapping (colliding) meshes, we enrich the notation of Logg [42] by the concepts 
of collision relations and collision maps. The collision relation To <->• 972 between 7o 
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Fig. 3.2. (Left) A typical example of a geometric predicate is the incircle test, which tells 
whether a given point d lies on the inside, on the outside, or on a circle defined by three points 
a, b, c. (Right) Axis aligned bounding boxes (AABB) of two triangles. 



and 972 is defined by 



% &T 2 = {(*, j) ■ Ti n F j + A T t € T Q A F 3 e dT 2 }. 



(3.6) 



The collision relation lists pairs of indices of all intersecting cells of the background 
mesh % with boundary facets of the overlapping mesh T 2 - Furthermore, each index 
of an intersected entity is mapped to the set of indices of intersecting entities via a 
pair of collision maps: 



We note that the collision maps To — > dT 2 and 972 — > To can be computed from the 
collision relation To <H- dT 2 - 

A naive approach to computing the collision relation between two meshes To and 
72 would be to intersect each cell of 7o with each cell of T 2 , resulting in an O(|7o|- l^l) - 
complexity, which is not feasible for large meshes. However, efficient algorithms and 
data structures which reduce the complexity dramatically have been developed in 
the fields of computer science, computational geometry, and computer graphics. The 
task of determining whether two objects collide arises naturally in the rendering of 
a computer game and the objects in question are often represented by meshes. In 
this paper, we limit the description to those techniques used in our work and refer to 
the books Ericson [27], Schneider and Ebcrly [57], and Akcnine-Moller et al. [3] for a 
broader overview. To efficiently find all pairs of the collision relation To dT 2l two 
important concepts are used: (i) fast intersection tests for pairs of simple geometric 
objects and (ii) spatial data structures for large sets of objects to accelerate collision 
queries. 

3.2.1. Fast intersection detection. Fast intersection detection for simple ge- 
ometric entities such as triangles and tetrahedra is based on so-called geometric predi- 
cates. Geometric predicates are tests which determine whether two geometric entities 
do intersect (collide) without computing the actual intersection. See Figure 3.2 for an 
example. This saves work since the actual intersection does not need to be computed 
when the geometric predicate is false. 



(To -+ dT 2 )(i) - {j : e To dT 2 }, 
(dT 2 -> 75)0') - : (i,i) e % » dT 2 }. 



(3.7) 
(3.8) 



Fig. 3.3. Bounding volume tree of gear mesh. (Left) Root bounding box (blue) and the next 
two generations (red and green). (Right) Every fourth level of the complete bounding box tree. 



3.2.2. Spatial data structures. Spatial data structures provide a hierarchical 
ordering of geometric objects, which allows quick traversal of large sets of geometric 
objects as part of a collision test. There exist two principal approaches based on cither 
space subdivisions or bounding volume hierarchies, which can both be realized as tree- 
like data structures. As the name indicates, subdivision approaches rely on some sort 
of geometric subdivision of the entire space embedding the structures of interest, in 
our case a finite element mesh. Typical examples are binary space partitions (BSP) 
trees, quadtrees, and octrees. Since the embedding space is subdivided, the leaf of 
a tree will usually not represent a single mesh entity and therefore further selection 
procedures are required to determine the actual intersections. 

In contrast, a bounding volume hierarchy is a tree which is built from bounding 
volumes; that is, simple geometric shapes containing the objects to be tested. Exam- 
ples of such data structures are axis aligned bounding boxes (AABB, see Figures 3.2 
and 3.3), oriented bounding boxes (OBB) and so-called k-DOPs, discrete orientation 
polytopes described by k hyperplanes. The purpose of using bounding volumes is 
twofold: (i) testing for collision of bounding volumes is cheaper than testing for col- 
lision of the bounded objects, and (ii) the simple geometry of the bounding volumes 
means that they can be stored efficiently in a hierarchical manner. 

3.2.3. Building the collision map. Bounding volume trees accelerate asym- 
metric collision queries between a tree embedding a large set of objects and a single 
simple object like a tetrahedron. The concept becomes even more powerful when 
intersection tests between two large sets of geometric primitives are desired, as in our 
case of two meshes. Then a hierarchical traversal of both trees greatly improves the 
C(|7ol ■ |72|)-complexity of the naive approach. Algorithm 2 describes how to employ 
a hierarchical traversal to compute the collision map. 

After the completion of the traversal according to Algorithm 2, which identifies 
all cells of 7o that intersect the boundary of 72, it remains to check whether the 
remaining cells (which do not intersect the boundary of I2) are either completely 
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Algorithm 2 Traversal of two bounding volume trees. Starting from the root of 
both trees, pairs of nodes are recursively tested for intersection. Only if the bounding 
volumes A and B of two nodes overlap may the children possibly overlap. A so- 
called descend rule is applied to decide with which of the two nodes to proceed. 
A popular and effective rule is to choose the one with the larger volume [27] since 
it gives the largest volume reduction for subsequent bounding volumes. If two leaf 
bounding volumes are reached, then the mesh entities bounded by them are tested 

for intersection. 

compute_collisions( J 4, B): 
if AC\B = 0: 
return 

else if is_leaf(^4) A isdcaf(B): 
if Ta H Tg 7^ 0: 

= (index(T4), index^)) 
% «■ dT 2 := To <-> dT 2 U 
else if descend_a(A, B): 
for a G children (A): 

compute_collisions(a, B) 

else: 

for b E childrcn(£>): 

compute_coilisions(A, b) 

descend_a(A, B): 

return isieaf(B) V (-.is_leaf(A) A \A\ > \B\) 



overlapped or not overlapped by T- To check this, it is enough to take a single point 
contained in each cell and check whether it is in f2 2 . To avoid building a third A ABB 
tree for the entire mesh Ti and to take advantage of the already built tree for 972, 
one can use a method known as ray-shooting; see Akenine-Moller et al. [3]. The idea 
is simple: counting how often a ray starting at the point in question intersects the 
surface dTi tells whether it is inside (odd number of intersections) or outside (even 
number of intersections). To efficiently find all ray-surface intersections, one may 
reuse the same AABB tree for 972 that was used in Algorithm 2. 

To summarize, the collision map and classification of cells according to the split- 
ting (3.3) can be found efficiently by utilizing AABB trees, fast traversal of pairs of 
AABB hierarchies, and ray shooting techniques. 

3.3. Mesh intersection. The next step is to compute the cut cells (polyhc- 
dra) {P; } and the interface decomposition {Tki}; see Figure 3.1. This computation 
may be phrased in terms of so-called boolean operations which are widely used in 
CAD systems to build complex geometries by performing boolean operations between 
primitives from a finite set of geometries. Algorithm 3 summarizes the use of boolean 
operations to compute {P ; } and {r fc; }. The resulting geometric objects are depicted 
in Figure 3.4. 

The boolean operations are completely delegated to the computational geometry 
library GTS [2] which uses a so-called ear-clipping algorithm to compute the surface 
tesselation of the intersection objects. The implemented ear clipping algorithms (first 
introduced in Fournier and Montuno [28]) are known to have 0(nlog(n)) complexity, 
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Algorithm 3 Mesh intersection. Note that the surface S defined by the union of all 
facets F intersecting the tetrahedron T divides T into two parts Tg and T<T according 

to orientation of S. 

for each T e To intersected by T: 
I = index(T) 

p i ■= T 
S:=$ 

for each F e (T -> dT 2 )(l): 

p? ■■= n 

for each F e 972: 
j = index(F) 
k = indcx(T F ) 
for eachTe {8% -> 

I = indcx(T) 

T k i :=fnr 



where n is the number of vertices of the 3D-polygon. This number in turn (n) depends 
on the mesh quality parameters such as the smallest and widest element angles and 
will in practice be bounded by a constant. 

3.4. Integration on complex polyhedra. The last required step in imple- 
menting the overlapping mesh method is to compute the integrals (3.4) and (3.5) 
on the cut cells {P/ 1 } and on the interface decomposition {Tki}, respectively. A 
widely used technique in the implementation of the extended finite element method 
[15, 24, 52, 58] is to decompose a cell into subcells which are aligned with the in- 
terface. The interface is either approximated by a plane within the cell [29, 30] or 
completely recovered [24, 58], even in the case of higher order interfaces [49]. Inside 
each tetrahedron of the subtetrahedralization, one may then use a standard integra- 
tion scheme. However, subtetrahedralization of an arbitrary polyhedron is in general 
a quite challenging problem [56] , and its existence cannot even be guaranteed without 
adding additional vertices [21], in contrast to the two-dimensional case. 

Methods for integration over arbitrary polygonal domains without the use of a 
subtriangulation have recently been presented in Natarajan et al. [53] and incor- 
porated into the extended finite element method with discontinuous enrichments in 
Natarajan ct al. [54]. Using Schwarz-Christoffcl conformal mappings as the funda- 
mental tool, the technique is strongly bounded to two space dimensions and hard to 
generalize to three space dimensions. 

We here propose an alternative approach, which is based on a boundary repre- 
sentation of the integrals. Despite the fact that this technique has been known for a 
long time, see for instance [40], it seems to be largely unknown within the finite ele- 
ment community. Our implementation is based on the efficient realization described 
in Mirtich [51]. This technique may be easily generalized to the computation of the 
general moment integral I a (P) over a polyhedron P, defined by 

I a {P)= j x a dx = J x^---x a d d dx, (3.9) 

where a = (ai, . . . , ay) denotes a multi-index of length d. 
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Fig. 3.4. (Left top) A tetrahedron cell T ; from the background mesh 7o overlapped by a tetra- 
hedral mesh Ti of a cube. (Bottom left) The resulting cut cell (polyhedron) P ; °. (Top right) An 
interface facet F S 972 of the overlapping mesh 75 intersected by a number of tetrahedra from the 
background meshTo- (Bottom left) The resulting facet decomposition (triangulation) {Tj.;}. 

The integral I a (P) is computed in three steps. The first step is to interpret the 
integrand as the divergence of a polynomial vector field and to rewrite (3.9) as surface 
integral: 

/ x a dx = / V-V-^ re l dx= V Vnf-e, / -f ds. (3.10) 

Jp Jp + p ^ p {r[ l J F d(a i + l) K ' 

Here, e, = (0, . . . , 0, 1, 0, . . . , 0) denotes the ith unit vector. 

Secondly, the plane equation ax + by + cz = d for each facet F allows the construc- 
tion of a projection map Z — h(X,Y), where X, Y, Z is some positive, orientation- 
preserving permutation of x, y, z. Using the parametrization h, one may rewrite the 
integrals of type f F x* ds in (3.10) as integrals in the XF-plane: 

fa^ds = --*-, [ x(X,Y,h(X,Y)fdXdY. (3.11) 

JF \ n Z\ Jh-^F) 

Here, nz denotes the Z-component of the normal vector rip. 

The third step consists of using a parametrization of dh~ 1 (F) and Green's the- 
orem in the plane to rewrite (3.11) as a sum of line integrals. Together, the three 
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steps reduce the evaluation of the moment integral (3.9) to the evaluation of a set 
of one-dimensional integrals with polynomial integrands. We note that the algorithm 
for the calculation of each moment integral has only a 0(#F) complexity. 

The moment integrals can be used in several ways to finally compute the inte- 
grals (3.4) and (3.5). First, one may compute the volume and the barycenter as the 
zeroth and first moments, respectively, which immediately provides a quadrature rule 
for exact integration of polynomials of degree one. Unfortunately, higher order mo- 
ments may not be directly reinterpreted as quadrature rules. Furthermore, as the 
construction of quadrature rules for higher degree polynomials can be quite involved 
[23], it can be advantageous to avoid run-time generation of quadrature rules. We 
therefore propose to use the moment integrals I a (P) directly by first interpolating 
the integrand onto a monomial basis and then summing the contributions: 



Ill-conditioning may be avoided for higher-order expansions by replacing the monomial 
basis by Bernstein polynomials. In the present work, linear Lagrange elements have 
been used throughout and so simple barycenter quadrature suffices. 

4. Implementation and data structures. We now discuss some of the data 
structures and classes which reflect the abstract concepts and algorithms described 
in Section 3. For some of these algorithms, we rely on existing implementations as 
part of the computational geometry libraries CGAL [1] and GTS [2], while other 
algorithms have been realized in the finite clement library DOLFIN [43, 46] as part of 
this work. Specialized algorithms for the Nitsche overlapping mesh method have been 
implemented as part of the extension library DOLFIN-OLM which is built on top of 
DOLFIN. The code is free/open-source, licensed under the LGPLv3, and available at 
http : //launchpad . net/dolf in-olm. 

4.1. The finite element library DOLFIN. Our implementation of Nitsche's 
method on overlapping meshes is based on the finite element library DOLFIN which 
is part of the FEniCS project [41, 44] for automated scientific computing. The main 
feature of FEniCS is the automated treatment of finite element variational problems, 
based on automated generation of highly efficient CH — h code from abstract high-level 
descriptions of finite element variational problems expressed in near-mathematical 
notation [4, 39]. This is combined with built-in tools for working with efficient rep- 
resentations of computational meshes [42] and wrappers for high-performance linear 
algebra libraries like PETSc [8-10] and Trilinos [36]. 

For this work, we have integrated the computational geometry libraries CGAL [1] 
and GTS [2] with DOLFIN. While CGAL and GTS provide a large part of the func- 
tionality needed to compute intersections between tetrahedra, the integration scheme 
was realized through an adapted version of Mirtich's code [51]. Furthermore, the 
assembly routine of DOLFIN was extended to handle integration over cut cells and 
meshes. Currently, the code for computing local interface integrals has to be imple- 
mented manually by the user, but we plan to extend FEniCS, in particular DOLFIN 
and the form compiler FFC [39, 45], to provide a full automation of Nitsche's method 
where a user only needs to supply the abstract variational problem (2.8). 

In the remaining subsections, we present the class abstractions of the algorithms 
and data structures described in Section 3. These new classes introduced in DOLFIN- 




(3.12) 
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void NitscheAssembler : : assemble_interf ace (GenericTensor& A, 

const NitscheFormfe a_nit , 
UFC& ufc_0, 
UFC& ufc.l) 

{ 

for (CutFacetPartlterator cut_f acet (overlapping_meshes) ; 

! cut_f acet . end() ; ++cut_facet) 

{ 

// Update quadrature cache to current facet part and return quadrature rule 
const QuadratureRule& quadrature = 

a_nit . interf ace_domain_quadrature_cache [cut_f acet->index()] ; 

// Get cells incident with this part of the cut cell 

std: :pair<const Cell, const Cell> cells = cut_f acet->adjacent_cells() ; 
const Cellfe cellO = cells. first; 
const Cellfe celll = cells . second; 

// Update ufc forms and update interface local dimensions 
uint local_facetl = celll . index(cut_f acet->entire_f acet ()) ; 
uf c_0. update (cellO) ; 
uf c_l .update (celll) ; 

// Tabulate dofs for each dimension on interface facet part 
// Tabulate interface tensor. 

a_nit . tabulate_interf ace_tensor (a_nit . interf ace_ A . get () , 

uf c_0 . cell , 
uf c_l . cell , 
local_f acetl , 
quadrature . size () , 
quadrature . points ( ) , 
quadrature . weights () ) ; 

// Insert matrix 

A . add(a_nit . interf ace_A . get () , interf ace_dofs) ; 

} 

> 



Fig. 4.1. C+ + implementation of the assembly of the so-called interface tensor accounting for 
the coupling between the two meshes 7i and 7~2- A CutFacetPartlterator provides iteration over 
all interface facet parts, represented by CutFacetPart and stemming from the facet decomposition 
{r^;}. Each CutFacetPart is associated with exactly one cell in each mesh, which can be accessed 
via the adjacent _cells member function. The overall design stresses the similarities to the assembly 
of facet contributions in the standard DC method. 



OLM allow the implementation of the Nitsche assembly algorithm in a descriptive and 
concise manner, as illustrated in Figure 4.1. 

4.2. The class AABBTree. The search data structure AABBTree has been added 
to the DOLFIN library. The implementation is based on the computational geometry 
library CGAL [1]. Basic search queries such as finding one or all cells intersecting a 
given entity or distance computation are exposed via an IntersectionOperator class. 
DOLFIN-OLM complements that functionality with providing a GTS-based AABB 
tree [2] to allow traversal of two bounding box trees as described in Algorithm 2. 
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4.3. The class OverlappingMeshes. The class OverlappingMeshes, provided 
as part of DOLFIN-OLM, is a key component in our realization of the overlapping 
mesh method. It mainly computes additional topological and geometric information 
to describe the overlap of the two meshes 7o and T2 and provides access to this 
information, in particular the collision relation 7o <-> 972 and the collision maps 
7o — > 972 and 972 — > 97o described in Section 3.2. 

A MeshFunction, as introduced in Logg [42], describes the splitting (3.3) of the 
mesh 7o by assigning different integer values to cells in the not overlapped part 7o,i, 
the completely overlapped part 7o,2, and the partially overlapped part 7o,r, respec- 
tively. In the same way, the boundary facets of the overlapping mesh 72 can be marked 
if the overlapping domain fi 2 is not completely contained in background domain f2o- 

4.4. Mesh iterators for overlapping meshes. Mesh iterators have been ad- 
vocated by Berti [16, 17] and Logg [42] and used among others by Bastian et al. [12] 
and Botsch et al. [18] as an important abstraction concept in mesh implementations 
and FEM frameworks [12, 43]. The iterator concept allows to access and iterate over 
mesh entities such as vertices, edges, facets and cells without knowing the details 
of the underlying mesh implementation. We have followed the same ideas in the 
case of overlapping meshes. The two pairs of classes CutCell, CutCelllterator and 
CutFacetPart, CutFacetPartlterator provide an interface to the intersected cells in 
7o,r and the interface facet partition -{Tj,;}, respectively. The FacetPart class and the 
corresponding iterator class mimic the original interface in DOLFIN by giving access 
to the two adjacent cells in the overlapping and the overlapped meshes, respectively. 
The interface assembly in Figure 4.1 presents a important use case. 

Similar iterator concepts have been used in Bastian et al. [13] where a general 
infrastructure to couple grids interfaced by the DUNE grid framework is presented. 
The CutFacetPartlterator corresponds to the Remotelntersectionlterator de- 
scribed in Bastian et al. [13] specialized to our case of Nitsche's method on overlapping 
meshes. Moreover, the CutCelllterator introduced here can be interpreted as an 
instance of Domainlntersectionlterator in Bastian et al. [13]. 

4.5. The class Quadrature and QuadratureRuleCache. The DOLFIN class 
Quadrature is a lightweight base class which only computes and stores quadrature 
data such as quadrature points, weights and order for a given polyhedron at run- 
time. Since the integration order depends on the underlying finite element scheme, 
the actual computation of the points and weights is meant to be implemented in 
subclass constructors. The class BarycenterQuadrature is such an instance of a 
subclass which computes a quadrature rule of order 2 for a given polyhedron based 
on the algorithm outlined in Section 3.4. For each intersected cell or facet part, a 
QuadratureRule object can be stored in a QuadratureRuleCache instance to save 
geometry and quadrature rule recomputation if several integrations have to be per- 
formed on the same intersected entity, as is the case for the computation of the stiffness 
matrix and load vector in the Nitsche overlapping mesh method. 

4.6. Forms and assembly. The DOLFIN Form form class represents the math- 
ematical concept of a finite element variational form. This class has been extended, 
as part of DOLFIN-OLM, to reflect the domain decomposition character of the over- 
lapping mesh method. A so-called NitscheForm class holds the description of the 
variational problem on each part fii and f2 2 of the domain. The coupling between the 
two forms is accomplished through a member function tabulate_interf ace_tensor, 
which computes the local interface tensor corresponding to (3.5). This is in addition 
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to the standard tabulate_tensor functions defined in the UFC code generation in- 
terface [5] for cell integrals, exterior facet integrals, and interior facet integrals (cf. 
Algorithm 3.1). In addition, the NitscheForm class gives access to the two over- 
lapping meshes as well as a quadrature cache to avoid recomputation of quadrature 
rules. 

Degrees of freedom of the cells of 7o that are entirely covered by the overlapping 
mesh 72 are inactive; that is, those degrees of freedom do not determine the solution 
and can be assigned an arbitrary value. For practical reasons, inactive degrees of 
freedom are included in the linear system but are set to zero by inserting 'one' on the 
diagonal of the corresponding rows and 'zero' in the right-hand side vector. This is 
automatically handled by the DOLFIN function Matrix: : ident_zeros. 

5. Numerical examples. To demonstrate the efficiency of our realization of 
the overlapping mesh method, we study two test examples. The first example is 
the Poisson equation. Here, an analytical solution for a suitable source function 
allows to verify the implementation by means of convergence studies. In addition, 
timings for the computation of both a standard Pi finite element method and its 
corresponding Nitsche approximation are presented. A further breakdown of the 
timings in the Nitsche case gives a clear picture of the additional costs associated with 
the geometry-related computations and their effect on the overall computation time. 
From that perspective, comparing Nitsche with a simple piecewise linear, continuous 
finite element method represents the most challenging test case, since then the extra 
work required for geometry-related computations contributes the most to the overall 
assembly time. 

As a second example, a linear elastic equation is solved with discontinuous mate- 
rial parameters at the interface between two overlapping meshes. 

5.1. Poisson equation. We consider the elliptic model problem (2.1)-(2.5) on 
the domain O = (0, l) 3 C K 3 . The overlapping domain £l 2 is a translation and 
rotation of the cube f2 2 = [0.3331, 0.6669] 3 according to Figure 5.1. The source 
function / is given by 

f(x, y, z) = 3(2tt) 2 sin(27ra;) sin(27n/) sin(27rz). (5.1) 

Homogeneous Dirichlet boundary conditions are used on the entire boundary. The 
exact solution is given by 

u(x, y, z) = sin(27rx) sin(27ry) sin(2-7rz). (5.2) 

The penalty parameter 7 in (2.9) is set to 7 = 50. A code extract from the imple- 
mentation of the solver based on DOLFIN-OLM is shown in Figure 5.1. 

The Nitsche approximation was computed on a sequence of meshes with decreas- 
ing mesh size h max — 1/N for the tessellations To and T2, starting at N — 14 and 
stopping at N = 104. For the integration over cut cells and interface facets, barycenter 
quadrature was employed. The resulting linear systems were solved by the precondi- 
tioned conjugate gradient method as implemented in PETSc in combination with the 
algebraic multigrid solver from Hypre used as a preconditioner. For both the stan- 
dard finite element method and the overlapping mesh method, a mesh independent 
number of CG iterations was observed (3-4 and 9-11, respectively). All computations 
were carried out on a Macbook Pro equipped with a 2.66 GHz Intel Core i7 processor 
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// Function spaces and forms on the overlapped domain 
Poisson3D_l: :FunctionSpace Vl(meshl); 
Poisson3D_l : :BilinearForm al(Vl, VI); 
Poisson3D_l : :LinearForm L1(V1); 
Ll.f = f; // assign source 

// Function spaces and forms on the overlapping domain 
Poisson3D_2: :FunctionSpace V2(mesh2) ; 
Poisson3D_2: :BilinearForm a2(V2, V2) ; 
Poisson3D_2: :LinearForm L2(V2) ; 
L2.f = f; // assign source 

// Build Nitsche forms 

PoissonNitsche : : BilinearForm a_nit(al, a2) ; 

PoissonNitsche : : LinearForm b_nit(Ll, L2, 

a_nit . overlapping_meshes_ptr () , 

a_nit . overlapped_domain_quadrature_cache_ptr () , 

a_nit . interf ace_domain_quadrature_cache_ptr () ) ; 

// Assemble 
Matrix A; 
Vector b; 

NitscheAssembler :: assemble (A, a_nit, 0, 0, 0); 
NitscheAssembler :: assemble (b, b_nit , 0, 0, 0); 

// Apply boundary conditions 
DirichletBC be (VI, uO, boundary); 
be . apply (A , b) ; 

// Solve linear system 
Vector x; 

solveCA, x, b, "eg", " amg_hypr e " ) ; 

// Split solution vector according to domains 
Function ul(Vl) ; 
Function u2(V2) ; 

a_nit .distribute_solution(x, ul, u2) ; 



Fig. 5.1. Code example for the Poisson problem. The domain decomposition character of 
Nitsche's method is clearly reflected by defining forms on each mesh separately and "gluing" them 
together via a NitscheForm. Since both forms are assembled into a single matrix, the solution has 
to be split via the distributesolution function. 



and 8 GB of RAM (1066 MHz DDR3). The benchmarks were repeated 10 times and 
averaged to obtain the reported results. 

Convergence. As the theoretical results recalled in (2.11)-(2.12) predict, an op- 
timal convergence rate in both in the H 1 - and L 2 -norm is observed; see Figure 5.4. 
Figure 5.3 clearly illustrates the smooth transition from the solution part u\ defined 
on the overlapped mesh 7i to the part U2 defined on the overlapping mesh T2 ■ 

Benchmarks. The comparison in Figure 5.5 shows that the Nitsche overlapping 
mesh method is only twice as expensive as the standard finite element method for 
the same mesh size and approximately the same number of degrees of freedom. A 
breakdown of the computing time shows that both the assembly and the solution of 
the linear system become twice as expensive when the overlapping mesh method is 

17 



Fig. 5.2. Mesh configuration for the Poisson problem. 




Fig. 5.3. Two-dimensional cross-section of the three-dimensional solution of the Poisson prob- 
lem showing good agreement between solutions computed on the overlapped and the overlapping mesh. 
(Left) Both solution parts patched together. The solution on the overlapping mesh can be seen far to 
the left. (Right) A zoom revealing the small discontinuity between solutions at the common interface. 

used. The latter can be attributed to the observed higher iteration numbers. The 
CPU time for the linear solve displays a kink in the slope at ca 3 and 4.5 million 
cells for the standard and Nitsche FEM, respectively, which may be attributed to the 
problem size increasing beyond a hardware-specific threshold (cache or memory size) 
that induces overhead. 

A further breakdown of the assembly time is shown at the bottom of Figure 5.5. 
We note that while the number of cells scales like 0(h~ 3 ) on a regular grid, the number 
of facets scales like 0(h~ 2 ). Consequently, all purely interface related operations with 
linear complexity should scale like iV 2//3 , where N is the number of cells. This is also 
the case for the interface decomposition, computation of the cut cells, and integration 
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Fig. 5.4. Convergence of the Nitsche approximation in the L 2 - and H 1 -norm for the Poisson 
problem. 

over the interface and the cut cells. On the other hand, the initial collision detection 
between the overlapped and the overlapping mesh is somewhat more expensive as it 
involves the searching of tree-like data structures (adding a logarithmic factor to the 
complexity). However, this extra cost is negligible compared to the total computing 
time. In conclusion, the timings depicted in Figure 5.5 indicate that the cost of 
assembly for the Nitsche overlapping mesh method is comparable to standard finite 
element assembly, and its relative efficiency increases with increasing mesh size. 

A note on ill-conditioned stiffness matrices. As analyzed in Burman [19] and illus- 
trated by numerical experiments in [47], the stiffness matrix stemming from Nitsche's 
method on overlapping meshes and related schemes can be ill-conditioned if the in- 
tersection gives very small elements, compared to the original clement size. Different 
approaches have been investigated to cure the schemes from ill-conditioned systems, 
either by choosing proper weights for the interface [38], replacing basis functions 
on small elements by extensions of basis functions on larger neighboring elements, 
or through the use of a so-called ghost-penalty [19, 20]. For the current work, the 
only precaution was to skip all cut cells with a relative measure of size smaller that 
10~ 15 . In our numerical experiments, we observed a mesh-independent iteration num- 
ber which indicates that ill-conditioning was not present. However, proper handling 
of small cells by introducing ghost penalties has been studied recently for the Stokes 
problem in Massing et al. [48]. 

5.2. Linear elasticity. As a second example, we consider a linear elastic body 
occupying a domain f^o in M 3 consisting of two subdomains f2j and Q2 with possibly 
different material parameters. The displacement and the stresses are assumed to 
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Fig. 5.5. (Top) Overall computing time spent by the standard FEM and the Nitsche FEM. (Bot- 
tom) Breakdown of the assembly time spent by the Nitsche method compared to the overall assembly 
time for the standard FEM. Note the expected slope of 1 for the standard (usual) assembly and of 
2/3 for the interface related operations. The grapWshows that, asymptotically, overall assembly time 
will be dominated by standard assembly. 



be continuous across the interface T = dtli fl c?S7 2 between the subdomains. The 
corresponding linear elasticity problem then takes the form: find the displacement 
u : fi -> such that 



div(<r(uj)) = L 


in f2j, i = 1, 2, 


(5.3) 


[<t(u) • n] = 


on r, 


(5.4) 


[u] = 


on r, 


(5.5) 


u = 


on 90o,d, 


(5.6) 


<r(u) • n = g 


on <9f2 0j N- 


(5.7) 



Here, the stress tensor a is related to the displacement vector u by Hooke's law 

c(uj) = 2/ijc(uj) + Ajtr(e(iij))I in Sl it i = 1,2, (5.8) 

where and /Xj are the Lame parameters in Oj for i = 1,2 and e(v) = (grad(v) + 
grad(v) T )/2 is the strain tensor. 

Nitsche's method for the Poisson equation (2.1) proposed by Hansbo et al. [33] 
can be adapted to the case of the linear elastic problem (5.3)-(5.7) and takes the 
following form: find u/j G such that 

a(u h ,v) = /(v) VveW,,, (5.9) 

where 

a(u,v) — \] / <r(ui) : grad(vi) dx (5.10) 

i=i 

-J a(ui)-n[v] dS- J <r(vi)-n[u] dS + ^ J h'^u] ■ [v] dS, 

^ v ' V ^ v ' S v ' 

Stress balance Symmetrization Penalty/Stabilization 

(5.11) 

l(v)= I f-vda;+ / g-vds, (5.12) 
Jn Jan 

with 7 a positive penalty parameter. Assuming that the elastic material is not nearly 
incompressible (A remains bounded), we may extend the analysis in Hansbo et al. 
[33] and prove optimal order a priori error estimates. See Hansbo and Larson [35] 
and Becker et al. [14] for details on discontinuous Galerkin methods for elasticity 
problems. 

Test configuration and numerical results. We consider the linear elasticity prob- 
lem (5.3)— (5.7) with fio = (—2, 2) 3 C M 3 which is overlapped by a propeller-like 
domain fti = (-1, 1) x P where P = (-1, 1) x (-0.2, 0.2) U (-0.2, 0.2) x (-1, 1) C M 2 ; 
see Figure 3.1. The right-hand side is zero and the boundary conditions on <9£1 are 
defined by 

'u = 0, on (-2,2) 2 x {-2}, 

<ff(u)-n = 0, on <9((-2,2) 2 ) x (-2,2), (5.13) 
ff(u) n = g on (-2,2) 2 x {2}, 
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Fig. 5.6. (Left) Undeformed, original domain consisting of a cube overlapped by a propeller-like 
domain. (Right) The deformed domain. The color bar corresponds to the norm of the displacement 
u. 



where g represents a combination of pure tangential, rotational force and a normal 
pressure: 

g(x,2/,z) = te|^~(0,0,2-y^MV) T . (5.14) 
5y x z + y z 

The Lame parameters are given by //, = Ei/(2 + 2z/j), At = 22; • Vi/((1 + Vi)(l — 2i^)) 
in Q.i for i = 1, 2, with E 1 = 10, E 2 = 0.1 • Ex, v 1= v 2 = 0.3. 

The numerical results are shown in Figures 5.6 and 5.7. The results indicate 
a "smooth" transition of the solution from the overlapped mesh to the overlapping 
mesh. 

6. Conclusions and outlook. We have demonstrated that overlapping mesh 
methods, in particular the Nitsche overlapping mesh method, may be implemented 
efficiently in three space dimensions through the use of tree search data structures and 
tools from computational geometry. Numerical tests show that optimal order conver- 
gence is obtained, that the overhead of the overlapping mesh method compared to a 
standard finite clement method is small (roughly factor two), and that the overhead 
is decreasing as the size of the mesh is increased. 

In the future, we plan to extend our implementation and the techniques studied in 
this work to handle fluid-structure interaction problems as well as contact problems. 
Furthermore, we plan to fully automate the implementation of Nitsche formulations 
on overlapping meshes by adding code generation capabilities for interface terms to 
FEniCS. 
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Fig. 5.7. (Left) Magnitude of the displacement u in a cross-section through the xy-plane of the 
solution from Figure 5.6. (Right) Magnitude of the displacement u in the cross-section through the 
xz-plane. Despite a very coarse resolution, both figures show a "smooth" transition of the solution 
from the domain Qi to Q2- 
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